DTIC ADA081328: Mathematical Model of a Mixed Surveillance System.

ADA 08

security clash ricAi-ioM or this rage r

REPORT DOCUMENTATIL

1. REPORT NUMlP

m

READ INSTRUCTIONS
FORE COMPLETING FORM

T I CATALOG NUMBER

i: title tmx&ifU) - * — - -

Mathematical Model of a Mixed Surveillance System

(. performing org. report number

a. performing organization name ano aooress

EPL Analysis

Olney, M3 20702 ^L/

'll. CONTROLLING OFFICE NAME AND AOORESS

CNO (OP-96) f/fj

Washington, D.C. 20350 '-''i

MONITORING AG§ilCNLNAM6_» AOCRESS<T/ dlHarant team Controlling Otliea) MS. SECURITY CLASS, (at thla raport)

0 J I UNCLASSIFIED

ISa. DECLASSIFICATION/ DOWNGRADING
SCHEDULE

16. DISTRIBUTION STATEMENT (at thlt R.pon;

Unlimited and approved for Public release.

19. KEY WOROS (Continue on rereroe ride if neceoomry and identify by block number)

Advanced Naval Vehicle Concepts Evaluation
ANVCE

Military Worth
Surveillance System

20. ABSTRACT (Continue on revert* tide If neceeeary m d Identify by block number)

This memorandum describes a model of reactive surveillance in which two
classes of contacts occur. A single search vehicle attsrpts to maintain seme
localization of a single target vehicle with contacts occurring inter¬
mittently as a Poisson process but only when the target is within a certain
range (exposure disk) of the searcher. An external surveillance system is
also present; it produces contacts as a Poisson process irrespective
of the target position. Searcher tactics are assumed to approximate optimal

DO ,2r» 1473 EDITION OF ( NOV 66 IS OBSOLETE

03

SECURITY CLASSIFICATION of this page (When DeteWntere*)

50

SECURITY CLASSIFICATION O F THU F AOEfRhan 0*

search based on an assumed cir yii nal target distribution

determined solely by the class " -Jfcher or external system)
of the most recent contact and < time since that contact.

SECURITY CLASSIFICATION OF This FAOEfWha*) Data Enfaratf)

£P£ jhJpi* ft *”**'<'

CONSULTANTS IN
applied mathematics, operations research

OWARO P. LOANS
TER S. SHOENPKLD
(301) •144B74

December 29, 1977

SUITS SIS

3*00 OLNCY-LAYTONSVILLS ROAD
OLNEY. MARYLAND SOS 3 2

MEMORANDUM

To:

ANVCE (OP- 96V)

Department of the Navy
Arlington, Virginia

Attn: Dr. H. L. Weiner

From:

Peter S. Shoenfeld

Subject :

Mathematical Model of a Mixed Surveillance System

i

This memorandum describes a model of reactive surveillance
in which two classes of contacts occur. A single search vehicle
attempts to maintain some localization of a single target vehicle
with contacts occurring intermittently as a Poisson process but
only when the target is within a certain range (exposure disk) of
the searcher. An external surveillance system is also present;
it produces contacts as a Poisson process irrespective of the
target position. Searcher tactics are assumed to approximate
optimal search based on an assumed circular normal target dis¬
tribution determined solely by the class (searcher or external

o£ the moSjt - recent

D'2VJ In

' ..nil D*

rf □

J U-: t : -

c r d [_

icat 1 o:i

■'7

~ 1 ~ i

f 1 1

'■Vi;: it.y Codes

Avail and/or

- 1 ‘it

k

special

.DISTRIBUTION STATEMENT F.

• , Approved lot public lelt'QM',
_ Distribution Unlimil< si

80 2 27 140

i

■t

■ V

The remainder of this memorandum is divided into sections.

The mathematical model and possible generalizations are
described in the first section. Numerical techniques employed
for the evaluation of limits and integrals are discussed in the
second. The computational algorithm is summarized in the third
section. A computer program listing is included as an appendix.

/V

» ' «

- <

Mathematical Model

Three types of contacts are considered; those by the
searcher, those by the external system when the target lies
within the searcher's exposure disk, and those by the external
system when the target lies outside this disk. Probabilistically,
the sequence of contact types and corresponding contact times is
described as a Markov renewal process in which the distribution
of next contact type and recontact time depends only on the type
of most recent contact; the process regenerates with each contact.
The sequence of contact types is a Markov chain. The intervals
between contacts are described separately as inhomogeneous,
continuous time parameter, Markov processes with different time
dependent transition functions and initial probability vectors
according to the type of most recent contact. These processes
have five states corresponding to the three types of contacts
(absorbing states) and target locations within and outside of the
searcher's exposure disk before recontact. The assumed circular
normal target distribution is described by a linear datum growth
law with initial datum area determined by contact type and datum
growth rate determined by target motion statistics.

These concepts are intended to model reactive surveillance
by ANVs; the external contacts might correspond to those generated
by SOSUS or by all other surveillance systems combined. However,
the concepts generalize to situations involving differing numbers
of contact types, transition functions, and datum growth laws.

The principal measure of effectiveness derived in this memo¬
randum is distribution of actual target range from the expected
(by the searcher) target location at a random time instant. This
has obvious application to anti-SLBM defense. Many other such
measures could be derived from the information generated by the
model.

- 3 -

The five states in the interval following a contact are
numbered and described as follows:

State 1 - Recontact by the external surveillance
system has occurred while the target is
exposed to contact by the searcher.

State 2 - Recontact by the external system has

occurred while the target is unexposed
to the searcher.

State 3 - Recontact by the searcher has occurred.

State 4 - Recontact has not yet occurred and the
target is exposed.

State 5 - Recontact has not yet occurred and the
target is unexposed.

The three types of contacts have the same numbering as the first
three states.

Let

1,2, 3, 4, 5
1,2,3

and

Pj(t) - (Plj(t),-**,P5j(t)), (j - 1,2,3)

where

t “ time since last contact.

Pij(t)

■ Probability

ity ^

process in
state i last
at time t contact of
after last type j ,
contact

Three separate Markov processes are obtained for j ■ 1,2,3. By
definition,

P^o) - (0,0,0, 1,0)

and

P2(0) - (0,0, 0,0,1).

Since the target must be exposed for contact by the searcher to
occur ,

P3(0) - (0,0, 0,1,0).

The model uses three absorbing states corresponding to
contact types although the use of only two (external and by
searcher) would appear more natural. This is done so that the
initial probability vectors, P-^(O), can be specified a priori.

In adapting the model to a more general situation in which there
are N (in this case two) classes of contacts and M (in this case
also two) intervening states, as many as MN (in this case three
sufficed) absorbing states might be required.

Let

5(A)

unconditional probability
that target is exposed to
searcher during search of
datum area A

and

*(A)

rate at which unexposed targets
enter exposure disk during
search of datum area A.

d

dT

Probability

target exposed
at time t + t

target
unexposed
at time t

1

T-0.

5

"Datum area" is defined by

A - 6tto2 (1)

2

where o is the parameter characterizing the assumed circular
normal target distribution. This is the area of a bounded region
which, if it contained the target with a uniform location distri¬
bution, would require the same expected effort to detect as is

2

required by a circular normal distribution with parameter o
(see reference Ca3). While in actuality the quantities § and »
depend on geographic details and are not functions of A alone,
the following formulas are reasonable approximations:

5(A) - 1 - exp(-n R£2/a) (2)

and

i(A) - 2ReVs/A,

where

(3)

R£ ■ radius of exposure disk, and
V - searcher patrol speed.

Formula (2) subsumes searcher allocation of effort and non¬
regularity in the search region. Formula (3) is familiar and
2

standard for ttR « a, it is shown below that it is sensible for
small A as well.

Consider an arbitrary Re and A. If the target is not
exposed to detection, it is confined to a region of expected area

- 6 -

Let “2 be the speed at which the boundary of the disc of exposure

sweeps through this unexposed area and let W be effective sweep
2

width. For fTRg « A, the entire boundary intersects A,

*S ■ V and W *■* 2Rg. For TTRg^ » A, only a small portion of the
boundary intersects A, and speed perpendicular to that segment
is V_( 2/n). But with probability .5, that segment is sweeping
into exposed area* i.e., no new area is exposed. Thus S ** S/tt
and W the length of boundary intersecting A.

Now if the target is located in the unexposed area* the
rate of entry into the disc of exposure is

SW

A exp (-ttRe /A)

Equating this to (3), one has

which, for ttRe « A, gives the expected result, W m 2Rg. As A
decreases in size* the right-hand member of (4) decreases* and
can be bounded by

W* (t5t)(£) (f^2 exp (-.5)) i

Taking *2 - S/n * a presumed minimum value,

(4)

7

(2.694),

W *

a very reasonable result. Thus formula (3) appears sensible
for small A, despite the inverse power dependence.

The assumed linear datum growth law after contacts of type

j is

Aj(t) " Aj(0) + pt (5)

where

Ai(t) ■ the datum area at a time t after a contact
J of type j, and

p ■ datum growth rate,

for j “ 1,2,3.

Since type 1 and type 2 contacts both come from the same external
surveillance system, the resulting initial datum is the same for
each (A^O) - A2(0)). The growth rate p is determined by target
motion statistics. If the target is moving in an unbiased random
walk with

- target speed, and

C “ mean time between target course changes,

it may be shown that the target location distribution is approxi¬
mately circularly normal with density in polar coordinates fr,0).

f(r,e)

2TTOZ(t)

exp

(6)

8

¥**-■

where

o 2(t ) - Vt2 C2 [t/C- (1- exp(-t/C))3.
Thus for large t,

So2(t) 9

9

and, since A * 6tto ,

p - 6TTVt2C . (7)

Formula (7) is used to compute p in applications.

Let

a - external recontact rate

d

&T

| Probability

external
recontact
time t + t

by

no recontact
by time t

»

T-0

X ■ searcher recontact rate for an exposed target

d

dr

Probability

searcher
recontact
t ime t + t

1#

by

target exposed
at time t and
no recontact
by time t

»

T-0

- 9 .

0i(t) ■ rate at which unexposed targets enter the
J exposure disk after a type j contact

Probability

target
exposed at
time t+t

target unexposed]
at time t and
last contact of
type j

*

T-0

and

Vj(t)

rate at which exposed targets leave exposure disk

< Probability

target
unexposed
t ime t + t

at

target exposed
at time t and
last contact of
type j

T=0

where

t ” time since last contact,

for j ® 1,2,3.

It is assumed that the rates a and X are actually constant
and that the rates pj and y j are functions of t alone. The model
requires these assumptions, which seem quite reasonable. It
follows that

dpj(t)

dt

PJ(t)

0

0

0

a

0

0

0

0

0

a

0 0 0

0 0 0

0 0 0

X -ta+\+Yj(t)] Y j (t)

0 0j(t) “[a>+p j(t)]

(8)

10 -

Now

Pj(t) - ’(Aj (t ) )

_ 2REVs
Aj(t)

W

by formula (3). For the case a - \ « 0, P^j(t) " §(Aj(t)) and
the fourth equation in (8) becomes

? ' (Aj(t) )Aj ' (t)

-Yj(t)S(Aj(t))+Pj(t)(l-5(Aj(t))).

Given a form for 5 this can be solved for Yj in terms of pj and
Aj. Using formula (2) for §,

Yj(t)

' exp(-HRE2/A1(t))
l-exp(-T'RE2Aj(t))

9j(t)

-

+

^ Aj (t) “J
at

CAJ(t)32

exp(-TTRE2/A.(t))

1 - exp(-TTRE Aj(t))

Pj(t) +

CAj(t)32

(10)

System (8) could be easily modified to discard external
contacts occurring at times when their information is inferior
to that already available from the last contact. This would
entail replacing the term a in (8) by

0 if A.(t) < A1(0)
a otherwise

This modification would have little effect in most cases.

11 -

Let

BtJ - Pjjtt) for i,J - 1,2,3.

It may be shown that since a > 0, the matrix tBjj3 is well
defined* has all positive entries, and all column sums one.

[Bij] is the transition matrix of the Markov chain characterizing
the sequence of contact types. It has a unique eigenvector
summing to one with eigenvalue one; this is the vector of steady
state probabilities. Let

<91>92’93> ’

steady state probability vector for
Markov chain characterizing sequence
of contact types

(ID

If

®13^®22 "" ^) “

23 12

(Bn - 1)(B22 - 1) - B12B21

Define

®23^®11 “ ^) “

13 21

(Bj^ ” 1)(B^^ - 1) - B-| <jBr

and

'22

12 21

*3 " ‘I

then

(12)

e, -

±

j *l+i2+i3

for j - 1,2,3.

Qj(t)

Probability

recontact within

recont

time t

last contact
of type j J

- 12

and

A j(t)

Qj(s)] ds .

Then

Qj(t) " i«! Pij(t) *

(13)

(14)

Integrating by parts,

Conditional expectation
of duration of portion
of interval between
contacts where time
since last contact < t
given that last contact
„of type j

I’* ( s if s s t'

J0 (t if s^t,

- fij(t).

(15)

In particular, defining

lim
t «* •

n:<c)-

(16)

Lj

Expected value of
time until recontact
after contact of
type j

(17)

This limit certainly exists since a > 0 and (8) imply that
l-Qj(t) is bounded above by exp(^xt).

13 -

for which the initial state is assumed known . Let

T . a random time instant, uniformly

[distributed on [0,T^] J *

i m the type of the last contact preceding
J1 T (0 if there is no such contact) J *

* m the type of the first contact following"] .

J2 T (0 if there is no such contact) J * ana

rp „ time elapsed at T since last contact
le (0 if there is no such contact) *

The random variables J2* and Tg all have limiting dis¬
tributions as which are independent of the initial state.

These distributions are:

Probability tjx - jJ - ^ (by^ll'.^n)! and (20))’

Probability [j2 - j] = iii (^({ij,

and (20);

(23)

Probability [Te < t]

S 8knk(t)

r \by (11), (15), and (20)/’

Now define

F(V) - Probability ?] , and

Aj'^Y) -

0 if Y s Aj (0)

unique t such that Aj(t) - Y otherwise

Here the model is regarded as a five state process with continu¬
ous time parameter in which the three states corresponding to
contacts are attained only for discrete instants.

15 -

Then by (24),

2 0^k(Ak’1(Y))
F(Y) - — -

r

With an assumed circular normal datum with parameter o , formula
(6) implies

target within
Probability distance R of
datum center

1 - exp(-R2/2o2)

Defining

Unconditional steady state probability
G(R) * that target is within distance R of ,
datum center

G(R)

Jf* 0 dF(6Trx)

I [ 1 - exp(-Rz/2x)3 - dx , using (1),

0 dx

Integrating by parts and changing variables,

G(R) - 3ttR

2 /** exp(-3TTR2/U)F(U)

Using formulas (25) and (5) for the functions F and Aj and
changing variables again leads to

G(R)

£ i Lf

r \^i\ kJ0

expC-3TTR2/(A^(0) +pt)3
~ (A^O+pt)*

nk(t)dt[. (26)

- 16 -

Evaluation of Limits and Integrals

System (8) is evaluated numerically by the Runge-Kutta
method for j ■ 1,2,3. The vectors P^(t) - (P^ j (t) , * * * , j(t))
and dP-Vdt - ' (t) , * * * ,Pg j ' (t)) are available explicitly at

each increment. The evaluation is carried out over the interval
C0* Tmax]» where Tmax is chosen so that the recontact probability,
QjCTjj^), will be quite large (say* .95) for each j. The contact
type transition probabilities are

Bij - 11--

Since

3

S B.. - 1
i-1 1J

necessarily, a convenient and reasonable approximation is

Bu~ VW* a1J’(Tmix) (x- J,pu<w) • (27)

The functions 0 j C t ) were defined as

flj(t) “ f [1-Qj(s)3ds
•'0
3

where Qj(s) - Pjj(s). Since values of the Pjj and their
derivatives are explicitly available at fine increments from the
Runge-Kutta evaluation of system (8), the 0^ can be obtained at
the same time to reasonable accuracy by trapezoidal integration

17

1

for times in the interval (0, . The limits

Lj - V0

are also required. Approximations to the Lj are obtained by
assuming that the probability of no recontact becomes a negative
exponential after a long time. This assumption is reasonable,
particularly if datum growth really ceases after a long time.
More precisely it is assumed that for t 2 Tfnax,

1-Q,(t)

- ] - * exp[-*.(t - T )3

1_ (T \ J tnax

max/

for some constant i y This leads to the approximation

x ^-Ql^max^2

L, ~ 0 ,(Tmov) + - 3 — — - ,

J j max 0 ' fT )

j ' max'

which is calculable since fij, Qj, and Q ^ ’ are all available for

t ~ Tmax*

The integrals in

G(R)

3TTPR"

I T k-1

tM

exp[-3TTR2/(Ak(0) + p (t)3

■2 - «k(t)dt

(Ak(0) +pt)'

(28)

must be evaluated for a number of values of R. The portion of
each integral in the interval (0, T_ov.] is evaluated from a table
of values of (t,fik(t)) which was saved from the Runge-Kutta
integration. The increment size changes several times in the

- 18 -

table* Simpson's rule is used separately on each Interval with
a constant increment size; trapezoidal integration is used
wherever an odd subinterval is left over. These integrals have
substantial "tails" in the interval [T •]. These are

IUclX.

estimated by assuming that for t > Tmax,

nk(t) ~ Lk, and
Ak(0) + pt - pt.

This leads to the approximation

3ttpIT 3

G(R) * - S l 0,

I T k-1

l

rmax exp[-37TR2/(Ak(0) +pt)]

(Ak(0)+pt)

2 - H^Ct) dt

+ l-expC-STTR^/pT^) .

(29)

19

Summary of Computational Algorithm

The algorithm is embodied in a computer program which is
included as Appendix A. Inputs are:

V. ■ searcher patrol speed (nm./hr.),

R£ “ radius of exposure disk (nm.),

a - external recontact rate (hr.”^) when target
is exposed,

X ■ searcher recontact rate (hr.”*0,

o

p * datum growth rate (nm./hr.),

Ai(0) ■ initial datum size after external contacts
(nm. ), and

A3(0) - initial datum size (nm. ) after searcher
recontacts.

The systems of linear first-order differential equations
(8) are Integrated over the interval [0, Tmny] by a library
Runge-Kutta routine for j ■ 1,2,3. The rates Pj(t) and Yj(t)
are calculated by formulas (9) and (10) using the datum growth
law (5). Values of Qj(t) and fij(t) are calculated simultaneously
as discussed earlier using formulas (13) and (14). Values of t
(time since last contact) and Oj(t) are retained for later use
for a number of increments. The reason for retaining the time
values is that the Runge-Kutta routine determines and adjusts its
own step size, defying user control. Values of t, Pjj(t),
p2j(t)» P3j(t), P5j(t), Qj(t), and Oj(t) are printed for

a number of values of t.

ft

The limits and Lj are calculated by formulas (27) and
(28). The steady state contact type probabilities are
calculated by formula (12). The Bjj, Lj, and 6j are then
printed out.

Values of the range distribution G(R) are calculated from
the tables of t, fij(t) an<* printed for the values R ** 5,10,

15, ’**,200 nm. , as discussed earlier using formula (29).

Peter S. Shoenfela

Reference [a]: L. D. Stone, Theory of Optimal Search. Academic
Press, New York, 1975.

21

APPENDIX A

COMPUTER PROGRAM LISTING

/

C

5

10

20

C

C

30

50

100

110

C

C

200

C

250

DIMENSION PRMT < 5 ) f DERP <5)fP<5)f AUX < 8 f 5 )

DIMENSION PH<3f3> fDPD<3f3> fT<3f400) f0MEGA<3f400>

DIMENSION OMINF < 3 ) f THETA (3)fB(3f3)fN(3)

REAL LAMDDA

COMMON VS » ALPHA f R r LAMBDA r RHO f AO 1 1 A02 f J » PD » DPD r T V OMEGA t
2 OMINF fNfTLAST fPDLAST
DATA PRMT/O. r 400. f 0.0625 f 0.0001 fO./

DATA NDIM/5/ fTMAX/300 • /

READ PARAMETERS
TYPE 10

FORMAT ( 1 X f ' VS t RE t ALPHA f LAMBDA f RHO f AZER01 f AZER03? ' / >

ACCEPT 20 f VS f R f ALPHA f LAMBDA f RHO f A01 f A02
FORMAT <7G>

EXTERNAL FCTfOUTP

INTEGRATE O.D.E. SYSTEM FOR EACH CONTACT TYPE
BY RUNGE KUTTA METHOD USING SUBROUTINE RKGS
PRMT < 2 ) =TMAX
DO 100 J=1f3
PRMT <5)=0.

DO 30 N=1 f 5
P ( K ) =0 .

DERP<K>=0.2

IF< J.NE.2)P(4>-1.

IF<J.EQ.2)P<5>=1.

N( J)=0
OMINF < J ) =0 •

TLAST =0 .

PDLAST =0 •

TYPE 50 f J

F0RMAT<20Xf' TYPE ' f 13 f ' CONTACT ' //4X f 'T' f SXf'PI'f 4Xf
2 ' P2 ' f 4X f ' P3 ' f 4 X f ' F" 4 ' r 4Xf'P5'f 4X f ' PD ' f 4X f ' OMEGA '// >

CALL RKGS ( PRMT f F* f DERP f ND I M f IMLF f FCT f OUTP f AUX )

TYPE IIOfIHLF

FORMAT ( IXf ' IHl.F=' f I4//>

FIND LIMITS TO OBTAIN TRANSITION MATRIX B(IfJ) AND
MEAN HOLDING TIMES OMINF<J)

DO 200 J=1 , 3

PDTOT =PD ( J f 1 ) +PD < J f 2 ) +PD < J f 3 )

DPDTOT =DF‘D ( J f 1 ) +DPD ( J f 2 ) +DPD < J f 3 >

OMINF < J)=OMINF< J)+( 1. -PDTOT >**2/DPDT0T
DO 200 1 = 1 f 3

B < I f J ) =PD < J f I ) +DPD ( J * I ) * ( 1 . -PDTOT > /DPDTOT

FIND EIGENVECTOR OF B(IfJ) TO GET STEADY STATE PROBABILITIES
D=(B<1f1)-1. >*<B(2f2>~1. )-B<1f2>*B(2f1)

THETA < 1)»(B<1f3)#<B(2f2)-1 . > -B < 2 f 3 > *B < 1 f 2 > > /D
THETA < 2 ) = <B<2f3)JK(B(1f1>-1.)-E«(1f3)#B<2f1)>/D
THETA<3)=-1 ♦

SUM=THETA< 1 )+THETA( 2 ) +THETA< 3 )

DO 2-0 J=1f3 r BRA-1

THETA ( J>=THETA< J)/SUM IHIS ^ lOUDC _ - ^

C PRINT THESE RESULTS

TYPE 260f0MINFf THETA* < <B(If J) » J-l»3> *1=1*3)

260 FORMAT ( //20X * ' MEiiAN HOLDING TIMES'//

2 3<3X*Er13.5> //20X * ' STEADY STATE PROBABILITIES V/3<3X»F13f5>//

3 20X * ' TRANSITION MATRIX ' //3< 3< 3X *F13 .5) / > )

C OAl.CtM ATI: RANOL HI SIR I HUT ION AND PRINT

TYI L 2/0

2/0 FORMA I < // 1 X * ' RANGE ' * SX r ' CUM . PROD . ' // )

DO 600 K«J1 * 40
RR»3»*K
PROD -*0.

D0 300 J«l»3
SUM-0.

NN=0

T0LD=0.

290 IF<NN.EQ.N( J) )G0 TO 400

H1=T ( J * NN+1 ) -TOLD

IF<NN+2.LE.N< J> >H2*T< J*NN+2>-T< J*NN+1>

IF < M2 . NE . HI . OR . NN+2 ♦ GT . N < J) )G0 TO 300

SUM=SUM+H1*<ELT (NN?RR)+4» KELT <NN+1 *RR > +ELT <NN+2*RR> >/3.

TOLD=T < J * NN+2 )

NN=NN+2
GO TO 280

300 SUM=SUM+H1*<ELT<NNfRR)+ELT<NN+1fRR> >/2.

TOLD-T < J * NN+1 )

NN=NN+1
GO TO 280

400 PR0B=PR0B+THETA< J)*SUM

500 CONTINUE

PROB-9 . 424778*RH0*RR*RR*PR0B/ < THETA < 1 > JKOMINF < 1 > +

2 THETA<2)*OMINF(2)+THETA(3>#OMINF(3) )+

3 1 . -EXP ( -9 . 424778*RR*RR/ < RHO*TMAX ) )

TYPE 510 »RR» PROD

510 FORMAT <1XfF6.1*2XfF7.3)

600 CONTINUE
GO TO 5
END

FUNCTION ELT < I »RR>

C EVALUATES INTEGRAND FOR RANGE DISTRIBUTION CALCULATION

DIMENSION PD < 3 * 3 > f DPD < 3 * 3 ) * T ( 3 * 400 ) t OMEGA ( 3 » 400 > »

2 OMINF ( 3 ) f THETA <3>fB(3f3)fN(3)

REAL LAMBDA

COMMON VS f ALPHA f R f LAMBDA f RHO f AO 1 f A02 f J f PD f DPD f T f OMEGA f
2 OMINFfNfTLASTfPDLAST
A0=A01

IF( J.EQ.3)A0=A02

IF<I.GT.0» AND «I.LE.N(J))GO TO 20
ELT ==0 •

GO TO 30
20 TT=T(JfI)

ELT=EXP<-9.424778*RR*RR/<A0+RH0*TT> >*OMEGA( JfI)/

2 <A0+RH0*TT>**2
30 RETURN

END

oi r j

I -

i

C

C

10

0

0

70

*

SUBROUTINE FCKTT. P.DERP)

EVALUATES DERIVATIVES FOR RUNGE-NUTTA ROUTINE
DIMENSION P<5>,DERP<5>

DIMENSION PD< 3 » 3 ) » DPD ( 3 » 3 ) t T ( 3 » 400 > t OMEGA < 3 1 400 > t
2 OMINF < 3 > » THETA ( 3)»B(3»3)»N<3)

REAL LAMBDA

COMMON VS » ALPHA > F< » LAMBDA , RHO » AO 1 » A02 » J , PD * DPD r T » OMEGA »
2 OMINF »N»TLAST »PDLAST
A(X)=AO+RHO*X
BETA< X )=2 . #R#VS/A ( X )

PP(X)=EXP<-3.141593*R*R/A<X> )

GAMMA<X)=<PP(X)/(1 .-PP<X> ) >*<BETA<X)+3.141593*R*R*RH0/A<X>**2>
A0=A01

IF (J.EQ.3)A0=A02
DERP<1)=ALPHA*P<4>

DERP < 2 ) ==ALPHA*P < 5 )

DERP < 3 ) =LAMBDA*P ( A )

DERP < 4 ) =- ( ALPHA+LAMBDA-f GAMMA ( TT > > *P ( 4 ) + BETA ( TT ) *P ( 5 )

DERP < 5 ) =GAMMA < TT ) *P ( 4 ) - ( ALPHA+BETA ( TT ) ) *P < 5 >

RETURN

END

SUBROUTINE OUTP< TT » P , DERP » IHLF » NDIM * PRMT )

OUTPUT ROUTINE FOR RUNGE-KUTTA INTEGRATION
REAL LAMBDA

DIMENSION P ( 5 ) » DERP ( 5 ) » PRMT < 5 )

. DIMENSION PD<3»3) »DPD(3»3) »T (3»400) »QMEGA(3r400> t
2 OMINF (3) >THETA(3) »B(3»3) »N<3>

COMMON VS > ALPHA » R » LAMBDA » RHO t A01 r A02 » J» PD » DPD r T r OMEGA r
2 OMINF »N»TLAST rPDLAST
, DO 10 1=1 , 3
F'D < J* I )=P< I )

DPD < J » I ) =DERP < I )

PDD=P < 1 ) +P < 2 > +P < 3 )

OMINF < J ) =OMINF < J ) + < TT-TLAST ) *< 1 . -(PDD+PDLAST) /2 . >

TLAST =TT
PDLAST=PDD

IF<AMOD<TT, ,25). GE. 0.0001 )G0 TO 50
IF < AMOD ( TT » 2 . ) . GE ♦ 0 . 0001 . AND . TT , GT . 4 . >G0 TO 50
IF<AM0D<TT»8.).GE. 0.0001. AND. TT.GT. 20. >G0 TO 50
IF ( AMOD < TT f 16 . ) • GE. 0 • 0001 • AND . TT . GT , 100 , ) GO TO 50
TYPE 20 f TT r P » PDD f OMINF < J >

FORMAT <lXrF7.2»6<lX»F5.3)flX»F6.2)

IF(AMOD(TTf ,0625) .GE. 0.0001 )G0 TO 70

IF(AMOD<TTr 1. ) .GE. 0.0001. AND. TT.GT. 4. >G0 TO 70

N ( J )=N < J ) +1

T(J?N(J>) =TT

OMEGA ( J »N < J > ) ®OMINF ( J )

IF (PDD.GT . 1 . )PRMT (5) = 1 .

RETURN

END lr

A- 3

lmniiOUl INI:. KKOil

c

c

c

o

i;

c

<:

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

r*

w

c

c

c

c

t - 1 IK I - 1 1 ; I

in um.vi;: a system ok first order ordinary differential
i;i.iiwiiimi with o i vi in inihae values.

USAGE

CALL RKOS (PRMT » Y .DERY .NDIM. IHLF.FCT .OUTP.AUX)
PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT.

DESCRIPTION OF PARAMETERS

PRMT - AN INPUT AND OUTPUT VECTOR WITH DIMENSION GREATER

OR EQUAL TO 5 » WHICH SPECIFIES THE PARAMETERS OF
THE INTERVAL AND OF ACCURACY AND WHICH SERVES FOR
COMMUNICATION BETWEEN OUTPUT SUBROUTINE (FURNISHED

BY THE USER) AND SUBROUTINE RKGS • EXCEPT PRMT(5)

THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE
RKGS AND THEY ARE

PRMT ( 1 ) - LOWER BOUND OF THE INTERVAL ( INPUT).

PRMT ( 2 ) - UPPER BOUND OF THE INTERVAL (INPUT).

PRMT ( 3 ) - INITIAL INCREMENT OF THE INDEPENDENT VARIABLE
(INPUT) .

PRMT ( A ) - UPPER ERROR BOUND (INPUT). IF ABSOLUTE ERROR IS
GREATER THAN PRMT (A). INCREMENT GETS HALVED.

IF INCREMENT IS LESS THAN PRMT<3> AND ABSOLUTE
ERROR LESS THAN PRMT ( A ) /50 . INCREMENT GETS DOUBLED.
THE USER MAY CHANGE PRMT ( A > BY MEANS OF HIS
OUTPUT SUBROUTINE.

PRMT ( 5 ) -

Y

DERY

NDIM

IHLF

NO INPUT PARAMETER. SUBROUTINE RKGS INITIALIZES
PRMT (3) =0. IF THE USER WANTS TO TERMINATE
SUBROUTINE RKGS AT ANY OUTPUT POINT. HE HAS TO
CHANGE PRMT ( 5 ) TO NON-ZERO BY MEANS OF SUBROUTINE
OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE
FEASIBLE IF ITS DIMENSION IS DEFINED GREATER
THAN 5. HOWEVER SUBROUTINE RKGS DOES NOT REQUIRE
AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL
FOR HANDING RESULT VALUES TO THE MAIN PROGRAM
(CALLING RKGS) WHICH ARE OBTAINED BY SPECIAL
MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP.
INPUT VECTOR OF INITIAL VALUES. (DESTROYED)

LATERON Y IS THE RESULTING VECTOR OF DEPENDENT
VARIABLES COMPUTED AT INTERMEDIATE POINTS X.

INPUT VECTOR OF ERROR WEIGHTS. (DESTROYED)

THE SUM OF ITS COMPONENTS MUST BE EQUAL TO 1.
LATERON DERY IS THE VECTOR OF DERIVATIVES. WHICH
BELONG TO FUNCTION VALUES Y AT A POINT X.

AN .INPUT VALUE. WHICH SPECIFIES THE NUMBER OF
EQUATIONS IN THE SYSTEM,

AN OUTPUT VALUE. WHICH SPECIFIES THE NUMBER OF
BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS
GREATER THAN 10. SUBROUTINE RKGS RETURNS WITH
ERROR MESSAGE IHLF^U INTO MAIN PROGRAM. ERROR
MESSAGE IHLF- 12 OR IHLF-13 APPEARS IN CASE
PRMT ( 3 ) -0 OR IN CASE SIGN ( PRMT ( 3 )). NE . SIGN (PRMT<2>-

PRMT ( 1 ) ) RESPECTIVELY.

- A-4 -

nooooooonononnoononnnonnonoorinnrjoooonooooonnoo

FCT - THE NAME OF AN EXTERNAL SUBROUTINE USED. THIS

SUBROUTINE COMPUTES THE RIGHT HAND SIDES DERY OF
THE SYSTEM TO GIVEN VALUES X AND Y. ITS PARAMETER
LIST MUST BE X * Y* DERY . SUBROUTINE FCT SHOULD
NOT DESTROY X AND Y.

OUTP - THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED.

ITS PARAMETER LIST MUST BE X»Y *DERY * IHLF »NDIM»PRMT •
NONE OF THESE PARAMETERS (EXCEPT * IF NECESSARY »

PRMT ( 4 > * PRMT ( 5 ) * . . . ) SHOULD BE CHANGED BY
SUBROUTINE OUTP. IF PRMT<5> IS CHANGED TO NON-ZERO»
SUBROUTINE RKGS IS TERMINATED.

AUX - AN AUXILIARY STORAGE ARRAY WITH 8 ROWS AND NDIM
COLUMNS.

REMARKS

THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM* IF

(1) MORE THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE
NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE
IHLF*=1 1 > *

(2) INITIAL INCREMENT IS EQUAL TO 0 OR HAS WRONG SIGN
(ERROR MESSAGES IHLF=12 OR IHLF=13>»

(3) THE WHOLE INTEGRATION INTERVAL IS WORKED THROUGH*

(A) SUBROUTINE OUTP HAS CHANGED PRMT ( 5 ) TO NON-ZERO.

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
THE EXTERNAL SUBROUTINES FCT ( X * Y * DERY ) AND

OUTP( X* Y, DERY * IHLF » NDIM* PRMT) MUST BE FURNISHED BY THE USER.

METHOD

EVALUATION IS DONE BY MEANS OF FOURTH ORDER RUNGE-KUTTA
FORMULAE IN THE MODIFICATION DUE TO GILL. ACCURACY IS
TESTED COMPARING THE RESULTS OF THE PROCEDURE WITH SINGLE
AND DOUBLE INCREMENT. 1

SUBROUTINE RKGS AUTOMATICALLY ADJUSTS THE INCREMENT DURING
THE WHOLE COMPUTATION BY HALVING OR DOUBLING. IF MORE THAN
10 BISECTIONS OF THE INCREMENT ARE NECESSARY TO GET
SATISFACTORY ACCURACY* THE SUBROUTINE RETURNS WITH
ERROR MESSAGE 3 1 ILF ^ 1 1 INTO MAIN PROGRAM.

TO GET FULL FLEXIBILITY IN OUTPUT* AN OUTPUT SUBROUTINE
MUST BE FURNISHED BY THE USER.

FOR REFERENCE* SEE

RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL COMPUTERS*
WILEY* NEW YORK/LONDON* 1960 * PP. 110-120.

THIS PAC-S Ir

.‘."■TTY FRACIILABLB
- -

A- 5 -

c

c

SUBROUTINE RKGS < PRMT f Y r BERY f NDIH t IHLF » FCT f OUTP f AUX >

C

C

C

C

C

c

c

c

c

c

c

c

c

BI MENS I ON Y < 1 > f DER Y < 1 > f AUX < 8 f 1 > f A < 4 > f B < 4 > f C < 4 > f PRMT ( 1 >
DO 1 I-l »NDIM

1 AUX<8fI)=.06A666A7*DERY<I>

X=PRMT ( 1 )

XEND-PRMT < 2 >

H=PRMT ( 3 )

PRMT <5)-0.»

CALL FCT< X f Y f DERY )

ERROR TEST

IF ( H* ( XEND-X ) >38f37f2

PREPARATIONS FOR RUNGE-KUTTA METHOD

2 A<1>=.5

A < 2) =.2928932
A<3)=1. 707107
A<4)=.16A66A7
B < 1 ) =2 .

B(2) = l .

B<3)=1 ♦

B(4)=2.

C<1>=.5
C(2)=. 2928932
C<3)=1. 707107
C < 4 ) = . 5

PREPARATIONS OF FIRST RUNGE-KUTTA STEP
DO 3 1=1 fNDIM
AUX < 1 f I ) =Y < I >

AUX ( 2 » I > =DERY < I >

AUX(3»I)=0.

3 AUX ( 6 f I ) =0 .

IREC=0

H=H+H

IHLF=-1

ISTEP=0

IEND=0

START OF A RUNGE-KUTTA STEP

4 IF ( (X+H-XEND )#M)7f6f5

5 H=XEND-X

6 IEND=1

7

8
9

RECORDING OF INITIAL VALUES OF THIS STEP
CALL OUTP <XfYf DERY f IREC f NDIM f PRMT >

IF (PRMT (5) )40f8f40

I TEST =0 . .

ISTEP=ISTEP+1

rbA'

jO VD1®

t

* A-6 -

m

nn no .on noon

C START OF INNERMOST RUNGE-KUTTA LOOP
J=1

10 AJ=A< J>

B J=B < J >

CJ=C< J)

00 11 I-lfNDIM
R1=H#DERY<I>

R2=AJ*<R1-BJ*AUX<6»I>>

Y< I )=Y < I >+R2
R2=R2+R2+R2

11 AUX(6r I )~AUX<6» I )+R2-CJ#Rl
IF < J-4 ) 12 » 15 > 15

12 J=J+1

IF ( J-3) 13fl4»13

13 X=X+.5*H

14 CALL FCT(X»YfDERY)

GOTO 10

END OF INNERMOST RUNGE-KUTTA LOOP

TEST OF ACCURACY

15 IF<ITEST)16»1A»20

IN CASE ITEST =0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY

16 DO 17 I»1>NDIM

17 AUX<4» I )=Y< I )

ITEST=1

ISTEP=ISTEP+ISTEP-2

18 IHLF=IMLF+1
X=X-H
H=.5*H

DO 19 I»1 »NDIM
Y ( I ) =AUX ( 1 » I )

DERY ( I ) “AUX ( 2 f I >

19 AUX<&> I )=AUX(3» I )

GOTO 9

IN CASE ITEST-1 TESTING OF ACCURACY IS POSSIBLE

20 IM0D=ISTEP/2

IF( ISTEP-IM0D-IM0D)21 »23»21

21 CALL FCT ( X > Y » DERY )

DO 22 I=1>NDIM
AUX <5»I)=Y<I)

22 AUX ( 7 » I > =DERY < I )

GOTO 9

COMPUTATION OF TEST VALUE BELT *

23 DELT=0.

DO 24 1=1 »NDIM

24 DELT=DELT+AUX ( 8» I ) #ABS< AUX< 4 » I )-Y < I ) >

IF (DELT-PRMT <4> )28»28r25

A- 7 -

non o o no

C ERROR IS TOO GREAT

25 IF< IHLF- 10 > 26 » 36 *36

26 DO 27 1=1 *NDIM

27 AUX (4*1) =AUX < 5 » I >

ISTEP=ISTEP+IST£P-4
X=X-H
IEND=0
GOTO 18

RESULT VALUES ARE GOOD

28 CALL FCT <X»Y »DERY>

DO 29 1 = 1 *NDIM
AUX < 1 r I ) -Y ( I )

AUX ( 2 » I ) = DERY < I )

AUX(3f I )=AUX(6f I )

Y ( I ) =AUX ( 5 » I >

29 DERY<I)=AUX(7f I)

CALL OUTP<X-Hr YrDERY * IHLF *NDIM*PRMT >
IF (PRMT < 5) >40*30*40

30 DO 31 I=1*NDIM
Y ( I ) =AUX (1*1)

31 DERY < I )=AUX<2» I )

IREC=IHLF
IF( I END >32 *32 *39

INCREMENT GETS DOUBLED

32 IHLF=IHLF-1
I STEP- I STEP/2

H

IF ( IHLF >4*33*33

33 IM0D=ISTEP/2
IF < I STEP- 1 MOD*- 1 MOD >4*34*4

34 IF ( DELT- . 02*PRMT < 4 > ) 35 * 35 * 4

35 IHLF=IHl.F-l
ISTEP=ISTEP/2
H=H+H
GOTO 4

RETURNS TO CALLING PROGRAM

36 IHLF=1 1

CALL FCT <X» Y»DERY)

GOTO 39

37 IHLF-12
GOTO 39

38 IHLF=13

39 CALL OUTP<X* Y*DERY » IHLF »NDIM*PRMT>

40 RETURN
END

Vh-

A-8 -